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1 INTRODUCTION 



X 



In the companion paper (|Ogilvie fc Latter] |2013| . hereafter 
Paper I), we have revisited the basic theory of warped astro- 
physical discs, in which the orbital plane varies with radius. 
We have introduced a local model, a warped shearing sheet, 
that is suitable for detailed analytical and numerical studies 
of warped discs. As we ll as rederivin g the nonlinear theory of 
laminar viscous discs (|Ogilvielll999l ) by a simpler route, we 
have shown how the local model can be used more generally 
to compute the internal torque that drives the large-scale 
evolution of the shape and mass distribution of the disc, 
thereby connecting the local and global dynamics. 

An important feature of a warped disc is that the orbit- 
ing fluid experiences a horizontal pressure gradient, antisym- 
metric about the midplane, that oscillates at the orbital fre- 
quency (or nearly so, if the shape of the disc is slowly varying 
in a non- rotating frame of reference). This force generates 
horizontal shearing motions, which are resonantly amplified 
if the disc is Keplerian (or nearly so), owing to the coinci- 
dence of the orbital and epicyclic frequencies. By transport- 
ing angular momentum in the radial direction, these motions 
provide an internal torque that causes an anomalously rapid 
evolution of the disc. 

Existing treatments of warped discs represent 
these internal motio n s, if at all, as laminar flows. In 
|Papaloiz"ou" fc Pringld (|l983l ) their amplitude is limited 
only by viscous damping, with the result that the internal 
torque, and the rate of evolution of the warp, are inversely 
proportional to the viscosity. (In this context, viscosity 
is usually taken to represent the effects of unresolved 
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physical processes such as small-scale turbulence.) In a 
propagating bending wave (|Papaloizou fc LinI Il995l ) the 
time-dependence ensures that the internal motions have a 
finite amplitude, even in the absence of viscosity, but the 
linear theory often predicts hypersonic flows above and 
below th e midplane. T he same is true in the nonlinear 
theory of lOgilvid (|l999l l. which we have revisited in Paper I. 

It has previously been noted that these 
oscillatory shear flo ws are likely to be lin- 
early unstable ('Papaloizou fc TerquemI Il995l : 

iGammie. Goodman fc O gilvie 2000). To the best of 
our knowledge, instability has not been detected in any of 
the global numerical simulations of warped discs (many 
of which are cited in Paper I), probably because they 
have insufficient spatial resolution or too much viscosity, 
either explicit or numerical, to allow the unstable modes to 
develop. Therefore these simulations also represent warped 
discs as laminar flows and may be missing an important 
physical process. 

Hydrodynamic instability is likely to lead to turbulence, 
or at least significant wave activity, and to alter the internal 
fiows in a warped disc. We expect an important modification 
of the large-scale dynamics of warped discs, and this is one 
reason why we are motivated to make a detailed study of the 
instability. Another reason is that it will provide a source 
of hydrodynamic activity in Keplerian discs, especially in 
poorly ionized regions where the magnetorotational insta- 
bility (MRI) is inoperative. These secondary motions could 
be important in determining many of the physical properties 
of these regions, for example by stirring up the dead zones 
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of protoplanetary discs. Note that only a very small warp 
may be required to excite these motions, owing to the reso- 
nant amplification of the internal flows in a Keplerian disc. 
Such small but dynamically significant warps may not be 
directly observable even with high-resolution imaging. More 
generally, we are interested in the interplay between hydro- 
dynamic instability resulting from a warp and other sources 
of turbulence and transport in astrophysical discs, including 
the MRl. 

The likelihood of a local linear hydrodynamic instabil- 
ity in warped discs was noted by IPapaloizou fc Tergueml 
and the first explicit calculations were 



made 



Gammie. Goodman fc Ogilvid 
l|2000l ) 



2003)EI 

iGammie. Goodman fc Ogilvie 1 20001 ') did not represent 
the warp itself, but considered the horizontal shearing 
motions, oscillating at the local orbital frequency, that are 
driven by the warp. These motions are linearly unstable 
when the associated shear rate is sufficiently large compared 
to the viscosity of the fluid, the instability taking the form of 
a parametric resonance of inertial waves. In nonlinear local 
simulations, the shearing motions decay as they lose energy 
to inertial waves. iTorkelsson et al.l (|2000l) studied the decay 
of similar flows in the presence of magnetohydrodynamic 
(MHD) turbulence driven by the MRL 

An effect related to the parametric instability of w arped 
discs is th e mode-coupling process investigated by iKatd 
l|2004 l and iFerreira fc Ogilviel l|2008l ). This allows a warp 
in the inner part of a disc around a black hole to excite in- 
ertial waves that are trapped near the radius at which the 
epicyclic frequency is maximal, and which may be related 
to observed high-frequency quasi-periodic osciUations from 
accreting black holes. 

In the local hydrodynam i c and MH D simulations of 
iGammie. Goodman fc Ogilviel (I2OO0I ) and ITorkelsson et al] 
l|2000h . the horizontal shearing motions were imposed as ini- 
tial conditions and were allowed to decay through instability 
or interaction with turbulence. However, in order to reach 
a better understanding of the role of instabilities or turbu- 
lence in the dynamics of warped discs, a model is required 
in which the effect of the warp can be maintained and the 
dynamics can reach a quasi-steady state. Such a model is 
provided by the warped shearing sheet deflned in Paper I. 

To prepare the way for future nonlinear simulations in 
the warped shearing sheet, we tackle the linear problem in 
this paper, analysing the hydrodynamic stability of the lam- 
inar flows comp uted in Paper I. Our analysis i s close ly re- 
lated to that of IGammie. Goodman fc Ogilvij (|200 0'). but 
is more detailed and is based on the nonlinear laminar flow 
solutions computed in a warped shearing sheet. We compute 
the growth rates of the unstable modes and their variation 
with several parameters. In certain regimes the instability 
can be clearly identifled as a parametric resonance of inertial 
waves, in agreement with previous work. We also speculate 
on the nonlinear outcome and its consequences for warped 
discs, as well as making suggestions for further investigation. 



^ The problem considered by iKumar Sz Colemanl l ll993l) . al- 
though motivated by warped discs, is too simplistic to be directly 
appUcable. 



2 LOCAL MODEL AND LAMINAR FLOWS 

In this section we summarize the equations of Paper I on 
which our analysis is based. The local model of a warped 
disc introduced in that paper is based on a circular reference 
orbit of radius ro and angular velocity f^o, where the dimen- 
sionless rate of orbital shear is q — — dlnSl/dlnr, equal to 
I for Keplerian orbits. The initial construction is identical 
to the standard shearing sheet or box and leads to the hy- 
drodynamic equations 



Dux — 2Q,o'Uy — 2qQ,QX — dxh, 

Duy + 2QoUx — ~dyh, 

Duz = —^l^z — dzh, 

Y)h = —cl{dxUx + dyUy -f dzUz), 

where 

T) = dt + Uxdx + Uydy + Uzdz 



(1) 
(2) 
(3) 
(4) 

(5) 



is the Lagrangian derivative and u is the velocity. For sim- 
plicity, we consider an isothermal gas for which 



P = CsP, 



(6) 



where Cs = constant is the isothermal sound speed, and we 
make use of the pseudo-enthalpy 



h = Cgln p + constant. 



(7) 



To represent a warped disc we introduce the warped 
shearing coordinates 



t = t, 



(8) 
(9) 
(10) 

(11) 



y = y + iTx, 

z' = z + cos T X, 
where 

r = not' = Qot (12) 

is the orbital phase and \'ip\ is the dimensionless warp ampli- 
tude. These coordinates follow the warped orbital motion. 
Partial derivatives transform according to 

dt = d't + qfloxdy ~ \il)\Qo siriT xd'z, (13) 

dx ^ d'x + qr d'y + \tP\ COST d'z, (14) 

dy = d'y, (15) 

dz = d'z, (16) 

so that 

B = d't + Vxd'x + {vy + qTVx)d'y + {vz + \tp\cosT Vx)d'z, (17) 
where 

Vx = Ux, (18) 



Vy = Uy + qfloX, 

Vz = Uz — IV'l^o sin r x 



(19) 
(20) 



are the relative velocity components. 

The hydrodynamic equations are therefore transformed 

into 



Dvx — 2Q.oVy = —{d'x + qTd'y + Ii/jI cost 



(21) 
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T)vy + (2 - g)f7of^ = -d'yh, (22) 

Duz + |'i/'|f2o sinTVx = —Qqz' — d'^h, (23) 

D/i = — Cs[(9^ + gr 9^, + [tpl cosTd'^)vx + d'yVy + d'^v^]. (24) 

These equations are easily extended to include a dynamic 
shear viscosity /i = ap/Q,o and a dynamic bulk viscosity 
Atb = QbP/f^o, where a and Qb are constant dimensionless 
coefficients. We do not write out the viscous terms explicitly, 
however. 

The simplest solutions of equations H21[) "(|24 p . extended 
to include viscosity, are laminar flows of the form 



Vx{z',t') = u{t)<:Ioz', (25) 

Vy{z',t') — v{t)Qoz', (26) 

v4z',t') =w{T)noz', (27) 

h{z',t') = clf{T)-^Qlz'^giT), (28) 



where u, v, w, f and g are dimensionless 27r-periodic func- 
tions that satisfy the equations 

drU + (w + \ip\ cos T u)u ~ 2v = \tp\ cos T Q 

— (ob + I Of) 1^1 COS T g{w + \ip\ COS r u) 

— a(7[|?/>| sinr + (1 + cos^ (29) 

drV + (lo + IV"! cos r u)j7 + (2 — q)u 

= — gl?/)! COST + (1 + IV"!^ cos^ ''")f]i (30) 

A-rW + (ui + cos r it)ui + sin r u = p — 1 

— (ob + \a)g(w + IV)] cosru) 

— a(7[|V'|^ sin r cost + (1 + cos^ t)w\, (31) 

d^/ = — (™ + It/"! cos r u), (32) 

drP = — 2(to + |?/)| cos r it)(7, (33) 

where dr denotes the ordinary derivative A/At. Numerical 
solutions of these equations are presented in Paper L 

3 STABILITY OF LAMINAR FLOWS 
3.1 Simple expectations 

Before embarking on a detailed analysis of the stability of 
the laminar flows, we make some very simple estimates, 
omitting factors of order unity. 

We are usually interested in discs that are nearly Kep- 
lerian (|q — |j <C 1) and of low viscosity (a <C 1). In such 
discs the horizontal laminar flows are resonant and can reach 
large amplitudes, even if the warp amplitude is small. 

The magnitude of u (or v, which is similar) can be sim- 
ply estimated as follows. The horizontal motion is a lin- 
ear oscillator, whose natural frequency is the epicyclic fre- 
quency. It is forced at the orbital frequency by the horizontal 
pressure gradients associated with the warp, and is damped 
by viscosity. The problem is similar to a damped harmonic 
oscillator forced at, or close to, its natural frequency; the 
response of the oscillator is limited either by detuning or 



damping, whichever is larger. When detuning of the reso- 
nance due to non-Keplerian rotation is the limiting factor, 
u scales as IV'l/k — |l (valid when this quantity is <1). 
When viscous damping of the horizontal flows is the lim- 
iting factor, u scales as IV'I/q (valid when this quantity is 
< 1) . Examples of these scalings can be seen in flgs 6 and 7 
of Paper I. 

These oscillatory shear flows are unstable in the absence 
of viscosity, and the inviscid growth rate of the instability 
scales with the shear rate uVl (|Gammie. Goodman fc Qgilvid 
|2000! ). Since the unstable modes have a non-trivial vertical 
structure and a horizontal wavelength related to the scale- 
height of the disc, H, viscous damping quenches the insta- 
bility if v/H^ > «f2, i.e. if a > it. Therefore (omitting factors 
of order unity) instability is expected when I?/)! > |q -- | |a 
or , whichever is larger. These are only very rough esti- 
mates or scaling relations; detailed results are obtained in 
Section 13.61 below. 

3.2 Linearized equations 

To examine the stability of the laminar flows we linearize 
the basic equations ((2T])"(l24]) about the solution (f25|) - (l28l) . 
introducing infinitesimal perturbations of the form 

Re[(5t)(z', r) exp(iA;ia;' -I- ifcyj/')], (34) 

etc., where (fe^, k'y) is a real horizontal wavevector. The lin- 
earized equations admit solutions of this plane-wave form 
because the basic equations (|21|| - (|24|) from which they are 
derived do not contain x' or y' explicitly, and the basic state 
about which we are linearizing is also independent of x' 
and y' . Using the relations ((9)1 and (|10|) between the primed 
and unprimed horizontal coordinates, we find 

exp(ifc^a;' -I- ik'yy') = exp(ifca;a:: -|- iky-y), (35) 

where 

{kx, ky) = {k'x + qk'yT, k'y) (36) 

is a wavevector that may depend on time. If fcy ^ then fea, 
depends linearly on r and we have a shearing wave, which 
corresponds to a non-axisymmetric spiral disturbance in the 
global geometry. If fcj, = then k^ = fc^ is a constant: the 
wave is axisymmetric and does not shear. 

We also adopt units in which f2o = 1 and Cs = 1, mean- 
ing that the unit of length is the scale-height Cs/Slo of an 
unwarped disc and the unit of time is Q,'^^ ■ The linearized 
equations in the absence of viscosity then take the form 

"DSv^ + [Sv^ -I- |?/)| cos r 5v.j:)u — 2 5vy 

= —{\kx + \ip\cosT d'^)5h, (37) 

'DSvy + [5vz + \'4>\ COST Svx)v + (2 — q)Svx = —ikySh, (38) 
DSv^ + {Sv^ + \t}}\ cos r Svx)w + \^\ sin r Sv^ — —d'^Sh, (39) 

DSh — {5vz + \ip\ COST 5vx)gz' = —ikxSvx 

— COST d'^Svx ~ ikySvy — d'^Svz, (40) 

where 

V = dr + i{kxU + kyv)z [w -\-\ip\ cos t u)z' d'^. (41) 
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For an axisymmetric wave {ky — 0), kx is constant and the 
coefficients of tfie linearized equations fiave a periodic depen- 
dence on T, with period 2n. For a shearing wave {ky 7^ 0), 
the coefficients are non-periodic because k^ depends linearly 
on T. 

The associated equation for the energy of the perturba- 
tion is (assuming suitable boundary conditions in z) 

d, / p {^\Svf + ^\Shf) dz' 

= — Re / p \Svxf\ip\cosT u + 5v*Svy(—q + l^plcosT v) 
+Sv^5vz{\i)\ sinr + \tp\ cost ui + u) 



+SvySvzV + ISv^l w 



dz'. 



(42) 



The right-hand side represents the exchange of energy with 
the laminar fiow and the orbital shear through Reynolds 
stresses. This equation can be used to bound the growth 
rate of the perturbation, as discussed in Appendix [Al 



3.3 Projection on to orthogonal polynomials 



To solve equations (|37|l - (|4ip we use a Galerkin spectral 
method, projecting the equations on to a basis of Hermite 
polynomials 



71 = 0,1,2,... (43) 



as used by lOkazaki. Kato fc Fukuj (|l987[ ) and many other 
authors: 



Svx{z' , t) — ^ it„(r) He„(z'), 

00 

5vy{z',T) = ^i;„(r)He„(z'), 

n=0 
00 

Svz(z',t) = ^ TO,i(r) He„_i(2') 

n=l 
00 

5h{z',T) = hnir) He„(z'). 



(44) 
(45) 
(46) 
(47) 



We further define Un — v„ — h„ — for n < and w„ = 
for n < 1. The Hermite polynomials, familiar from the 
quantum harmonic oscillator, are particularly well suited to 
linear waves in isothermal discs, having a Gaussian weight 
function that is proportional to the density of an unwarped 
isothermal disc. Using the relations 



He^s) = n}ie„-i{x), 

x'lie„{x) = Re„+i{x) + nHe„_i(a;), 

we find 



(48) 
(49) 



drUn + i{kxU + kyV)[Un-l + {n + l)u„ + i] 

+ (io -I- COST u)[rwi„ -I- (71 -I- l)(n -f 2)u„+2] 
+u(w„+\ + cos rti„) — 2vn = —ikxhn 

COST (n + l)/i„+i, (50) 



drVn + i{kxU + kyv)[v„^i + {n+ l)Vn + l] 

+ {w + cosru)[nVn + {n + l)(n + 2)w„+2] 
+v{wn+i + \ip\ cosrun) + (2 - q)u„ = —ikyhn, (51) 

drW„ + i{kj:U + kyv){w„-i + nw„+i) 

+ {w + \ip\ COST M)[(n - l)w„ + n{n + l)w„+2] 

+ w{Wn + l-ipl COSTUn-l) + sinT M„„l = —uhn, (52) 

drhn + i{kxU + kyv)[hn-i + {n+ l)hn + l] 

+ (w + cosru)[nhn + {n + l){n + 2)h„+2] 
^(.9 ^ i){i>'n + (n + l)w„+2 

+ \ip\ cost[ii„_i + {n+ l)u„+i]} 

= —\kxUn + COSTM„_l — \kyVn -f Wn. (53) 

Equations (|50p . (|5ip and (|53|l are to be solved for n > 0, and 
equation (|52p for n > 1. The equivalent equations including 
viscosity are given in Appendix [B] 

3.4 Solutions in the absence of a warp 

In the absence of a warp, \'4!\ = u = v = w = g — 1 = Q and 
the inviscid equations reduce to 



dr^n 2Vji — ikxhn, 
drV„ + (2 — q)un = —ikyhn, 

drWn = ~nhn, 

drhn = ikxUn ikyVn -f" Wn, 



(54) 
(55) 
(56) 
(57) 



with no coupling between different values of n. The same 
property holds when viscosity is included, but only if ab ~ 
|a. For axisymmetric waves [ky = 0), solutions cx e~"^'^ 
exist, with angular frequency lu given by, in the inviscid case. 



i-Lj'^ + n) [-u^ + 2(2 - q)] - fc^cj^ = 0, (58) 

(cf. lOkazaki. Kato fc Fukudl 19871 ) and, in the viscous case, 
a more complicated quartic equation. For n > 1, the high- 
frequency solutions of equation H58p have an acoustic charac- 
ter, while the low- frequency solutions have an inertial char- 
acter. In the case n = 0, however, the high-frequency so- 
lution is the well known two-dimensional (acoustic-inertial) 
wave, while the low-frequency solution is a zero-frequency 
vortical mode. The dispersion relation is illustrated in Fig.[T] 



3.5 Solutions in the presence of a warp 

In the presence of a warp, the different Hermite modes are 
coupled by numerous terms in equations (|50p - (|53p . For ex- 
ample, the terms proportional to i{kxU + kyv), which derive 
from the shearing of the waves by the laminar flow, couple 
neighbouring modes in such a way that an infinite ladder 
is produced. Indeed, the linearized equations (|37p - (|40p have 
no solutions that are strictly polynomial in z' in the presence 
of a warp, because the i{kxU + kyv)z' term in C generates an 
infinite power series in z' . Another type of coupling is pro- 
vided by the terms proportional to {g — 1) in equation H53p . 
These terms arise because the scale-height of the warped disc 
is time-dependent and different from that of the unwarped 
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Figure 1. Dispersion relation for axisymmetric waves {ky = 0) in an unwarped disc in the absence of viscosity, for q = 1.5 (left) and 
q = 1.6 (right). The angular frequency ui of the waves is in units of Qq, and the radial wavenumber kx is in units of f2o/cs- Different 
colours indicate different values of the vertical mode number n (0: black, 1: red, 2: green, 3: blue, 4: grey, 5: pink, 6: orange, 7: gold, 
8: magenta, 9: purple, 10: maroon, 11: cyan, 12: navy); only modes with n < 9 are shown for q = 1.5. Vertical bars of length 1 show 
the possible resonant couplings between neighbouring modes. These occur at wavenumbers for which two waves exist with vertical mode 
numbers n and n + 1 and with angular frequencies u} and ui + 1. 



disc on which the Hermite polynomials are based. For com- 
putational purposes, we truncate the system at some order 
A'" by setting u„, v„, Wn and h„ to zero for n > N. 

The parametric instability, at least in its simplest form, 
involves a three-mode coupling between the warp (together 
with the associated laminar flow) a nd two inertial waves 
l|Gammie. Goodman fc Ogilviell2000l ). The theory is worked 
out in detail in Appendix [C] which is based on an expan- 
sion of the equations for small In the local model, the 
warp and its associated laminar flow have radial wavenum- 
ber 0, vertical mode number 1 and frequency 1. In order 
to couple to the warp, the inertial waves should therefore 
have the same radial wavenumber as each other, while their 
vertical mode numbers should differ by 1 (owing to equa- 
tion I49p and their frequencies should add up to (or differ 
by) 1. The density of the spectrum of inertial waves means 
that there are infinitely many wavenumbers at which this 
coupling is possible; several of these are illustrated in Fig. [T] 
(Couplings of acoustic and inertial waves are found not to 
lead to instability.) 

We investigate only axisymmetric waves (ky = 0) in 
this paper, looking for instabilities of the laminar flows. As 
noted above, in this case the linearized equations have coef- 
ficients that are 27r-periodic in r. Our truncated system is 
a set of 4:N — 1 ordinary differential equations (ODEs) and 
we analyse their solutions using Floquet theory. The char- 
acteristic solutions (eigenmodes) of the system of ODEs are 
27r-periodic functions of r multiplied by e"^, where s is a 
complex growth rate. Choosing a set of 4A'' — 1 linearly in- 
dependent initial conditions (by setting all variables except 
one to zero at r = 0) and integrating the ODEs from r = 
to r = 27r using a Runge-Kutta method with adaptive step- 
size, we evaluate the monodromy matrix that determines 



the evolution of an arbitrary initial condition over one pe- 
riod. We calculate the eigenvalues of the monodromy matrix, 
which are the characteristic multipliers e^'^" , and deduce the 
complex growth rates. Growing solutions are found when 
Re(s) > 0, and indicate linear instability of the laminar 
flow. (The imaginary part of s gives the angular frequency 
of the eigenmode, but this is only defined modulo 1.) 

This numerical method has been verified by reproducing 
the solutions of the algebraic dispersion relation ((58]) and its 
viscous equivalent in the absence of a warp. In the presence 
of a warp, as described below, it produces results that agree 
with the analytical calculations of parametric instability in 
Appendix [C] Equations (|50p - (|53|) were also solved with a 
pseudospectral method, which provides a separate numeri- 
cal check on the results. As noted above, the characteristic 
solutions are 27r-periodic functions of r multiplied by e^"^, 
and if this common factor is extracted from the solutions, a 
term involving the complex growth rate s appears in each 
equation. After truncating the number of Hermite modes 
to A'', the T domain is partitioned into M segments and 
the T derivatives approximated using Fourier collocation. 
This yields a square matrix of size AMN, the eigenvalues 
of which are the complex growth rates s. Good agreement 
was achieved with the monodromy approach described in 
the previous paragraph. 

3.6 Numerical results 

We begin a discussion of the numerical results by consider- 
ing the relatively simple case of an inviscid non-Keplerian 
disc with q — 1.6 and a small warp of amplitude <^ 1. 
We compute the growth rate of the fastest growing mode for 
a range of radial wavenumbers and plot the results in Fig. [2] 
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Instability occurs in bands of k^c centred on the wavenum- 
bers predicted by the analysis of three-wave coupling; the 
red points in the figure show the growth rate at the centre 
of each band (up to kx = 5) calculated in Appendix ICll 
These growth rates are proportional to and tend to a 
limit as fcx — >■ oo. The bands of instability have a non-zero 
width, also proportional to l-i/)], because the parametric in- 
stability permits some detuning. The first few bands are 
well separated, but eventually they overlap and merge, pro- 
ducing a continuous curve. The limiting growth rate in the 
continuous regime is twice the limiting value predicted by 
a single three-wave coupling, because each mode (n) can 
engage simultaneo usly in three-wave couplings with it s two 
neighbours (nil) (|Gammie. Goodman fc Ogilviel2000 l). For 
larger values oi \tp\ the theory developed in Appendix[Cl] be- 
comes inaccurate. 

In Fig. [3] we show an example of a growing eigenmode, 
within the first band of parametric instability in the top row 
of Fig. [2] This mode consists mainly of a superposition of 
the n — 1 and n = 2 inertial waves, coupled resonantly by 
the warp so as to cause exponential growth. The resulting 
motion is neither symmetric nor antisymmetric about the 
midplane, but involves the oblique motions typical of iner- 
tial waves below the cpicyclic frequency (also seen in fig. 1 
of iGammic . Goodman & Ogilvie 2000). The variation of the 
velocity field with orbital phase is due to the interference 
of the two inertial waves with different phase speeds u)/kx; 
the pattern is not significantly distorted by shear, because 
both the warp amplitude and that of the laminar fiow are 
quite small. By referring to equation (|42|) and to the first 
panel of fig. 6 of Paper I (which shows that u is approx- 
imately proportional to sinr) we can understand how the 
mode grows: Svx and 5vz are negatively correlated at phase 
r = Ti/2 (panel 2) where the shear rate sinr + u reaches 
its positive maximum, and are positively correlated at phase 
r = 37r/2 (panel 4) where the shear rate reaches its negative 
minimum. The Reynolds stresses associated with this mode 
are therefore able to release energy from both the warped 
orbital motion and the laminar fiow. 

The situation for a Keplerian disc is more complicated. 
Viscosity must be included into order to obtain laminar flow 
solutions, and the analysis of parametric instability is mod- 
ified, mainly because of the different phase relationships in 
the laminar flows. Fig. |4] shows some results for very small 
warps and Fig. [5] for larger, but still probably unobservable, 
warps. The viscosity parameter a is chosen in each case to 
be sufficiently large that the laminar flow has a modest am- 
plitude {\u\ < 1). Distinct bands of parametric instability 
are clearly visible only for very small and a, as in the 
left panel of Fig. |4l where good agreement is found with the 
analytical calculation in Appendix IC2I Note that viscous 
damping reduces the growth rate at larger values of k^. For 
larger warp amplitudes, the analytical theory becomes inac- 
curate and instability occurs in broader ranges of k^. 

In situations where the laminar flows have large am- 
plitude (|it|^l) and the viscosity is small (a <C 1), large 
growth rates (Re(s) > 1) may be expected. However, in this 
regime the eigenvalues returned by our numerical method 
do not converge as the truncation order A'' is increased, 
and therefore we do not present any such results in this pa- 
per. Initially the difficulty takes the form of spurious modes, 
presumably arising from the truncation of the Hermite ba- 



sis, with non-convergent growth rates that can exceed those 
of the genuine modes. Such unwanted solutions can be ex- 
cluded by comparing the eigenvalues calculated with differ- 
ent truncation orders. Later, however, the entire spectrum 
fails to converge. The reason for this behaviour is probably 
that the coupling terms proportional to ik^u in equations 
(|50|) - (|53p are strong, and indeed the problem is worse for 
larger values of k^- It may be simply that the Hermite basis, 
which is designed for wave modes in an unwarped disc, is in- 
capable of handling the modes of a warped disc with strong 
laminar flows. 

It is also possible, however, that modal solutions cease 
to exist when the laminar flows are too strong, or when the 
radial wavenumber is too large. After all, it is not obvious 
why an unwarped isothermal disc supports acoustic and in- 
ertial waves whose energy is conflned near the midplane. 
If the vertical gravitational acceleration were constant, in- 
stead of being proportional to z, there would be a continuous 
spectrum instead of conflned modes. Similarly, in a stably 
stratified vertical isothermal disc with adiabatic exponent 
7 > 1, the internal gravity waves are not confined but form 
a continuous spectrum. In a shearing sheet without radial 
structure or boundaries, non-axisymmetric eigenmodes do 
not exist because of the inexorable shear, and are replaced 
by non-modal solutions. It is possible that the oscillatory 
shear of a strong laminar flow in a warped disc may prevent 
the formation of axisymmetric eigenmodes. 

So far we have discussed the main bands of instability 
that occur for radial wavenumbers kx^l. There are, how- 
ever, weaker forms of instability at longer wavelengths. Close 
to kx = 0.925 in Fig. [21 for example, is a band of instabil- 
ity (most easily seen in the bottom left panel) whose height 
and width scale with \ip\'^ rather than \Tp\, indicative of a 
weaker, higher-order mode coupling. We also find viscous 
overstability at sma ller values of fcj,. Visco us overstability 
(|Katdll97 3; see also iLatter fc Ogilviell2006l ) tends to cause 
a slow growth of long 'density' waves in a viscous disc. In 
the absence of a warp, the n = mode has a growth rate 
scaling as akx for kx I, although this can be suppressed 
by the addition of a sufficiently large bulk viscosity. Modes 
with n > 1 are not unstable. Overstability appears in our 
calculations for viscous warped discs, although it is not vis- 
ible in Figs 2] and [S] because the growth rates are too small 
compared to the parametric instability. In the absence of a 
warp, overstability occurs ai kx <^ 1 in an isothermal disc 
when ctb < |a in the case q — 1.5, or when ah < fa in 
the case q — 1.6. This illustrates the fact that overstability 
is more likely closer to a black hole, where (effectively) q is 
closer to 2. Our calculations suggest that the presence of a 
warp and fast laminar flows acts against overstability. For 
example, when q = 1.5 and ctb = 0, overstability is present 
at = but is suppressed for IV"! > I.Iq. 

We summarize our numerical results in Figs [6] and [71 
where we plot contour lines of the maximum growth rate 
obtained for < fc^, < 5 in the parameter space of warp am- 
plitude (li'l) and viscosity (a), for non-Keplerian {q — 1.6) 
and Keplerian {q = 1.5) discs. In the non-Keplerian case 
(Fig. [HJ parametric instability is found for a < 0.2226|?/'| in 
the limit <C 1, as determined in Appendix ICll and indi- 
cated by the dotted straight line. For larger the onset of 
parametric instability deviates somewhat from this straight 
line. The weaker form of instability found in most of the 
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Figure 2. Linear stability of laminar flows for q = 1.6. The growth rate of the fastest growing mode (in units of Qo) is plotted versus 
the radial wavenumber (in units of Qq/cb). The red points show the growth rate at the centre of each band of parametric instability (up 
to kx =5), based on the analysis in Appendix [CT] Larger growth rates are obtained where the resonance bands overlap. The left panels 
require resolutions of up to A'^ = 30 for convergence. The right panels, which illustrate the behaviour for larger wavenumbers, require up 
to AT = 90. 
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Figure 3. Example of an inviscid growing eigenmode, within the first band of parametric instability in the top row of Fig. [2] Velocity 
vectors &v, multiplied by p^^"^ to show the distribution of wave energy, are shown in the x' z' plane at over two radial wavelengths, at 
orbital phases r = 0, 7r/2, tt and 37r/2. 



remainder of the diagram is the viscous overstabihty, which 
favours longer wavelengths. Fig. [7] for the Keplerian case 
covers a smaller range of warp amplitudes. Parametric in- 
stability is found for a < 0.123Q\tp\^^^ in the limit < 1, 
as determined in Appendix IC2I and indicated by the dot- 
ted parabola. Again, viscous overstabihty provides a weak 
growth of longer waves above this boundary. 



In the lower right corner of Fig. [3 the strength of the 
laminar flows means that growth rates larger than 0.1 are 
expected, but we have not computed this region systemati- 
cally because of the onset of a lack of numerical convergence 
as described above. For similar reasons we do not extend the 
contour plots to larger values of li/il and a. 



4 DISCUSSION 

To the best of our knowledge, the linear hydrodynamic in- 
stability of warped discs that we have described has not 
been seen in any of the global, three-dimensional numeri- 
cal simulations of warped discs. This is probably because 
the simulations have insufficient spatial resolution or too 
much viscosity (explicit or numerical) to allow instability 
to devel op. For example , the h igh-resolution SPH simula- 
tions of iLodato fc Pried (j2010l ). whic h appea r com patible 
with the laminar nonlinear theory of lQgilvi3 (|l999l ). have 
estimated shear viscosities in the range 0.04 ^ a ^ 0.5, in 
addition to bulk viscosity. Their largest (initial) values of 
\ip\ are 0.026 for their low-amplitude warp and 1.35 for their 
high-amplitude warp. Although we would expect instabil- 
ity to occur when the amplitude is high and the viscos- 
ity is low, it is not clear whether the unstable modes can 
be meaningfully resolved, or escape catastrophic damping. 
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Figure 4. Linear stability of laminar flows for q = 1.5, including viscosity. The left panel is computed with A'^ = 20 and the right panel 
with TV = 40, which are sufficient to provide converged results for this range of kx. The red points show the growth rate at the centre 
of each band of parametric instability, based on the analysis in Appendix IC2I and taking into account viscous damping; this analysis is 
inaccurate for the right panel. Larger growth rates are obtained where the resonance bands overlap. 
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Figure 5. Continuation of Fig. [5] to larger warp amplitudes. Af = 30 is sufficient for the left panel and N = 80 for the right panel. 



when the SPH smoothing length is typically 0.6H . Further- 
more, it is in this regime that the warp tends to steepen 
into a 'break'. The high- resolution grid-based simulations of 
iFragner fc NelsonI (|2010l ) have a grid spacing of 0.2H in the 
vertical direction and an explicit shear viscosity in the range 
0.005 < a < 0.1. The warps are mostly very small, however. 
The most warped disc they simulate is also the most viscous, 
and the instability is probably absent for that reason. 

Given the inherent limitations of current global simula- 
tions, it would make sense to study the nonlinear outcome of 
the hydrodynamic instability of warped discs with local nu- 
merical simulations in the warped shearing box as defined in 
Paper I. Such simulations should be able to determine the 



more complicated flows that replace the laminar solutions 
in the unstable regions of the parameter space, and allow a 
calculation of the torques acting in a warped disc in the pres- 
ence of the instability. Here we can only speculate on what 
the outcome might be. Since the instability feeds from the 
shear in the laminar flows, it is possible that it will reduce 
the amplitude of those flows and may therefore limit the as- 
sociated torques that determine the evolution of the warp, 
as discussed in section 2 of Paper I. In this way the rapid 
diffusion of warps in Keplerian discs might be ameliorated, 
allowing warps to survive more readily. 

Although we have considered only axisymmetric insta- 
bilities {ky = 0) in this paper, non-axisymmetric shearing 
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Figure 6. Contours of growth rate, maximized over < < 5, 
for q = 1.6 and = 0. Dashed contours: 0.0005, 0.001, 0.0015, 
0.002, 0.0025, showing the region of viscous overstability. Sohd 
contours: 0.01, 0.02, 0.03, . . . , 0.12, showing the region of para- 
metric instability. Dotted line: a = 0.2226\ip\, the analytical pre- 
diction for the onset of parametric instabihty for <C 1. 



waves are also likely to be amplified, if only transiently. The 
nonlinear evolution could be computed initially with two- 
dimensional simulations (independent of y), although the 
behaviour in three dimensions may be different. 

It will also be important to determine the interaction of 
a warp and its associated hydrodynamic instabihty with the 
magnetorotational instability (MRI) in discs that are suffi- 
ciently coupled to a magnetic field. This can be done by solv- 
ing the MHD equations in a warped shearing box. If the MRI 
is sufficiently powerful, it may suppress the hydrodynamic 
instability; the presence of a magnetic field will certainly 



modify the low- lrequency par t 
in the disc (e.g. lOgilvid [l99i ). 



However, if the warp and 
its associated flows are sufhciently strong, their effe cts are 
bound to dominate over the MRI. Indeed, Torkclsso n et al.l 
were able to make a numerical study of the interac- 
tion between the MRI and the horizontal shearing motions 
associated with a warp, and found evidence that the hy- 
drodynamic instability dominates if the shear is sufficiently 
strong. 

The model we have investigated in this paper is inten- 
tionally simplified by being isothermal. If the disc is stably 
stratified, buoyancy forces will modify the details of the in- 
stability by replacing the inertial waves with inertia- gravity 
waves. However, since buoyancy forces vanish at the mid- 
plane, their effect may be limited. 
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Figure 7. Contours of growth rate, maximized over < fc^: < 5, 
for q = 1.5 and 0^ = 0. Dashed contours: 0.00005, 0.0001, 
0.00015, . . . , 0.0003, showing the region of viscous overstability. 
Solid contours: 0.01, 0.02, 0.03, . . . , 0.1, showing the region of 
parametric instability. Dotted line: a = 0.1236|</'|^/'^, the analyt- 
ical prediction for the onset of parametric instability for \^\ <^ 1 . 
The region of fast laminar flows in the lower right part of the 
diagram is excluded from the calculation. 



5 CONCLUSION 

Even warps that are too small to be observable may have 
important dynamical consequences for astrophysical discs. 
The main reason for this is the coincidence of the orbital 
and epicyclic frequencies in a Keplerian disc. A warp that 
is stationary, or nearly so, in a non-rotating frame, implies 
oscillatory pressure gradients that drive horizontal shearing 
motions close to resonance. The simplest laminar fiow solu- 
tions provide a large torque that leads to the anomalously 
rapid diffusion or propagation of warps in Keplerian discs 
that has been noted in previous work. However, these flows 
are often linearly unstable and will be replaced by more com- 
plicated solutions and associated torques that remain to be 
computed. This dynamics will alter the evolution of warped 
discs. It also has the potential to introduce non-trivial hy- 
drodynamic behaviour in the form of turbulence, or at least 
wave activity, in discs with very small warps, especially if 
no other form of turbulence is present. This effect could be 
important, for example, in activating the dead zones of pro- 
toplanetary discs if there is any small misalignment in the 
system such as that caused by a planet or a distant binary 
companion on an inclined orbit. 

We have used a local model, the warped shearing sheet, 
to compute the laminar flows and to analyse their linear sta- 
bility in Keplerian and non-Keplerian warped discs. Hydro- 
dynamic instability derives from the parametric resonance of 
inertial waves and is widespread, especially when the viscos- 
ity is low or the rotation is c lose to Keplerian. Our analysis 
is clos ely related to that of iGammie. Goodman fc Ogilvid 
l|2000t ). but is more detailed and is based on the nonlinear 
laminar flow solutions computed in a warped shearing sheet. 

Future work on warped discs ought to include numerical 
simulations of the nonlinear outcome of the hydrodynamic 
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instability in the warped shearing box in two or three dimen- 
sions, and also of its interaction with the magnetorotational 
instability. The linear analysis presented in this paper should 
provide some initial guidance for this work. 

Closely related problems exist for eccentric discs and 
tidally distorted discs. In both cases the non-circular stream- 
lines present the fluid with an oscillating geometry that 
leads to the parametric excitation of inertial waves. These 
are versions of the elliptical i nstabil ity, well known in fluid 
dynamics (e.g. iKerswelUbOOd ). Thus lGoodniaiil (|l993l ) anal- 
ysed the elliptica l insta bility of tidally distorted discs and 
IRvu fc GoodmanI (| 19941 computed its nonlinear ev o lution 
in a two-dimensional local model, while iPapaloizoul (|2005l ') 
analysed the elliptical instability of eccentric discs. In each 
case the level of hydrodynamic activity that can be expected 
is limited by the ellipticity of the flow. An important differ- 
ence in the case of warped discs is that coincidence of the 
orbital and epicyclic frequencies in a Keplerian disc leads to 
a resonant enhancement of the internal flows, so that even 
a very small warp may produce significant hydrodynamic 
activity. 
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APPENDIX A: BOUND ON THE GROWTH RATE 

The right-hand side of equation (|42p can be written as the integral of an Hermitian form: 

~ J p{5vy ■E{T)Svdz', (Al) 
where E(r) is the (symmetric) rate-of-strain tensor with (dimensionless) components given by 

Exx = cos Tu, Eyy^O, E^^='w, (A2) 

Exy ^ Eyx = \{—q+ \ip\cosTv), Exz = Ezx = \{\il}\sinT +\'4!\cosTw + u), Eyz = Ezy = \v. (A3) 

Let the ordered eigenvalues of E(r) be E\{t) < E2{t) < Es{t). Then we have Ei\Sv\^ < (Sv)* -ESv < Es\Sv\'^ and so the 
quantity (|A1|) is bounded above by —Ei J p\5v\^ dz' . It follows from equation (|42|l that the instantaneous growth rate of the 
amplitude of the perturbation (in the energy norm) is limited by 



dr In 



1/2 



<-£i(r). (A4) 



For the axisymmetric waves discussed in Section [3.51 the Floquet growth rate is therefore bounded by the orbital average, 

s < -{El)r. (A5) 

In the absence of a warp the eigenvalues of E are ±g/2 and 0, leading to the familiar bound s < \q\/2, i.e. half the shear 
rate. In the presence of a warp, the eigenvalues do not have simple expressions, but various algebraic bounds can be placed 
on them if desired. An example of such a bound is 

Ef < EijEij = {\^\ cosTuf + ^{-q+ cosTvf + KlV"! sinr + cost w + uf + . (A6) 



APPENDIX B: VISCOUS LINEARIZED EQUATIONS 

When viscosity is included, the linearized equations (|37|l - (|4ip become 

TySvx + ((5uz + IV'I cos r 5vx)u — 25vy = — (ifc^, + IV"! cos t d'^)5h 

+ [ikx + \i>\ cosT{d'z - gz')][2a{ikx + cosTd'.^)5vx + (cib - fa) A'] 
+ (ifca; + IV'I cosr9^){[2Qj-i/jj COST It + (ab — |a)(™ + IV'I cost u)\5h} 
+ikya[{\kx + IV'I COST d'^)5vy + ikySvx + (— g + IV'I cosTv)5h] 

+a{d'z — gz')[{\kx + IV'I cost d'z)5v2 + d'^Svx] + a(|V'l sinr + IV'I cost w + ujO'^Sh, (Bl) 

"DSvy + {5vz + IV'I COST Svx)v + (2 — q)5vx = —ikySh 

+ [ika: + IV'I cost(92 - 5^')] [(i^^^ + IV'I cosTd'^)5vy + ikySvx]} 
+ {ikx + IV'I cos t di)[a{~q + \ cos r v)Sh] 
+iky{2aiky5vy + (ob — |a)[A' + [w + \tp\ cosTu)Sh]} 

+ {9'z — gz')[a{iky5vz + d'^Svy)] + avd'^Sh, (B2) 

VSvz + (Svz + IV-'I cos T Svx)w + IV'I sinr Sv^ = —d'^Sh 

+ [ikx + IV'I cosr(a!. - gz')]{a[{ikx + IV'I cost 9^)5^2 + d'^Sv^]} 

+ {ikx + IV'I cos r 9z)[q(|V'I sinr + IV'I cos tw + u)Sh] + ikya{iky5vz + d'^Svy + v5h) 

+ {d', - gz')[2ad'jv, + (ob - |a)A'] + [2aw + (ab - |q)(w + IV'I cosTu)]d'jh, (B3) 

VSh - {Sv, + IV'I cos T 5vx)gz = - A', (B4) 
with 

T> = Ot + i{kxU + kyv)z' + (ui + IV'I cos r u)z'd', (B5) 
and 

A' — {ikx + IV'I COST d2)5vx + ikySvy + d'^Sv,. (B6) 



© 0000 RAS, MNRAS 000, 000-000 



Hydrodynamic instability in warped discs 13 



When they are projected on to the basis of Hermite polynomials, they become 

dru„ + i{ka:U + kyv)[un-i + {u + l)u„+i] + {11) + l^p] COS T u)[nu„ + {u + l)(n + 2)u„+2] 

+ u{Wn+l + \lp\ COSrUn) — 2Vn = -ik^hn - jV"! COS T (jl + l)/ln + l 

+ikx{'2a[ikxU„ + COS r(n + + (ob - |ci)A„} 

cosT[g - l)(n + l){2a[\kxUn+-i + cos T(n + 2)u„+2] + {ctb - |a)A„+i} 

cosrg{2Q:[ifc2,u„_i + \-^\ cosrnu„] + (ob - |a)A„_i} 
+ [2Q:|'i/;| COST u + (ob — f + j^l cost u)][ikxhn + \4>\ cosr(n + 
+ikya[ikxV„ + \tp\ cos r{n + l)u„+i + ikyUn + {—q + cos r v)h„] 
-a{g — l){n + l)[ikxWn+2 + |V'l cos r(n + 2)to„+3 + {n + 2)un+2] 

~ag{\kxWn + \'4>\ COST nwn+i + nun) + ci{\tp\ sinr + \tp\ costw + u){n + l)/in+i, (B7) 

drVn + i{kxu + kyv)[vn-i + {u + l)vn+i\ + (w + ji/'l COS T u) [ni;„ + {n + l)(n + 2)v„+2\ 
+v(w„+\ + COS rii„) + (2 — q)un — —ikyh„ 
+ikxa[ikxVn + |^| cos r(n + l)v„+i + ikyUn] 

-\ip\ cos T{g - l)(n + l)a[ikxV„+i + cos r(n + 2)u„+2 + ikyU„+i] 
-\tp\ cos T ga[ikxV„-i + \tp\ cos Tnv„ + ikyU„^i] 
+a{—q + \ip\ cosTv)[ikxhn + lip] cosr(n + 
+iky{2aikyVn + (ob — |a)[A„ + {w + \ip\ cost u)hn]} 

— {g - l)(n + l)a[iky'Wn+2 + {n + 2)v„+2] — ga(ikyW„ + nv„) + av{n + l)h„+i, (B8) 
drWn + i{kxU + kyv)(wn-i + nWn+i) + (w + COST u)[(n — l)wn + n{n + l)w„+2] 

+ w{Wn + llpl COSTUn-l) + {ipl sinTUn-l = — 7l/l„ 

+ikxa[ikxW„ + \ip\ cosTUWn+i + nu„] 

-\^\ cosr(g - l)nQ:[ifc^iu„+i + ji/'l cosr(n + l)w„+2 + (n + l)u„+i] 

— 1^1 COST ga[ikxWn-i + \-tp\ cosT{n — + (n — 

+Q!(|'(/'I sinr + COST TO + u)[ikxhn-i + Itp] cosTnhn] 

+ikya{ikyWn + nvn + vhn-i) ~ (g — l)n[2a(n + l)w„+2 + (ob — |a)A„] 

—g[2a{n — l)wn + (ab - |a)A„_2] + [2aw + (ob - + IV'I cosTu)]nh„, (B9) 

drhn + i{kxu + kyv)[hn-i + {n + l)/i„+i] + {w + cosTu)[nh„ + (n + l)(n + 2)/i„+2] 

--g{w„ + {n+ l)u;„+2 + IV"! cosr[it,i_i + (n + = -A„, (BIO) 

with 

A„ = ikxUn + \ip\ cos r(n + l)u„+i + ikyVn + {n + l)wn+2- (Bll) 



APPENDIX C: PARAMETRIC INSTABILITY 
CI Non-Keplerian case 

Parametric instability of a slightly warped disc occurs when the warp p rovides a weak nonlinear coupling be tween wave modes 
of an unwarped disc that satisfy an appropriate resonance condition jGammie. Goodman fc Ogilvidlioool '). This is a special 
case of three-wave coupling in which one of the waves, namely the warp, is of much larger scale and is therefore effectively of 
zero wavenumber. A small amount of viscous damping of the waves can be allowed for, to compete with the growth provided 
by the parametric resonance. To analyse this regime we consider an expansion of the equations for small and small a, 
letting Q = \^p\ai and = |^|abi, where ai and abi are of order unity in the limit \tf}\ <C 1. To avoid the additional 
complications of the coincidence of the orbital and epicyclic frequencies in a Keplerian disc, we assume here that g 7^ §. The 
Keplerian case is discussed in Section [021 below. 
The laminar flow solutions have the expansion 

M= |V'|5'sinr + 0(|7/>|^), V ^\iP\CcosT + 0{\ij\'^), w ^ 0{\tpf), g - 1 = 0{\rljf), (CI) 

with 
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Then the hnearized solutions have the expansion 
«„(r) = u^°\ro, n, ■ • ■ ) + Mnl'\ro, n, ■ • ■ ) + O(l^l') 



(C3) 



and similarly for v„, Wn and /i„, where (to, ti, • • ■ ) = (r, IV'I''", ■ ■ ■ ) are multiple time-scales. Accordingly, the time-derivative 
At becomes do -\- \'4>\d\ -f ■ ■ • , where di — d/dri. The reason for allowing for multiple time-scales is that, in the regime of 
interest, the growth due to parametric resonance and the damping due to viscosity are weak and occur on a time-scale that 
is long, by a factor 0(|?/;|), compared to the orbital time-scale. In this analysis we are not interested in still longer time-scales 
and will suppress any dependence on T2, etc. 

We consider axisymmetric waves {ky = 0), which have a well defined frequency in the absence of a warp and are therefore 
suitable for parametric resonance. At the leading order, Odi/il"), the linearized equations yield 



LnUn ' 



0, 



where L„ and t/i"' are a linear operator and a vector of unknowns, given by 



do 


-2 





ihx 


2-9 


do 














do 


n 







-1 


do 



^ 71 



■ (0)- 

(0) 



(C4) 



(C5) 



These are just the linearized equations for axisymmetric waves in an unwarped disc. As discussed in Section 13.41 solutions 
exist that involve only a single value of n and are proportional to exp(— icjro), where the frequency lu satisfies the dispersion 
relation (|58|) . The corresponding eigenvector is 



U 



(0) 



nkxio 



n) 



At the next order, Od-i/)!^), the linearized equations yield 
r ttW _ 

where the effective forcing vector is 



-diU. 



(0) 



-ikxS sinr[u,*"2i + (n + - Ssmrw'^l^ ~ cosr(n + -|- Xn 

-ik^Ssinrlvll'l^ + {n + -Ccosrwl^'l^ + y„ 

-ikxS smT[w'^l-^ + nwl^li] - sin ru^^l^ + Zn 



'ikxS sinrlhl^l^ + ( 



,(0) 



and the viscous terms are 

Xn = -ai[{2kl + n)ui^^ + ik^wi^^] + (obi - |Qi)ifc,[ife,ui°^ + {n + l)™i°^2], 



Yn = -ai[{kl + 



„(0) 



+ qik^h(^>], 



~ai[{kl + 2(n - 1))™!"' - ifc.nu^f)] - (abi - faOlifc^ui^la + (n - l)wl^^]. 



(C6) 



(C7) 



(C8) 



(C9) 
(CIO) 
(Cll) 



In these equations, terms due to the warp couple each mode (n) with its neighbours (n ± 1). (The viscous terms also couple 
n with n ± 2 if ab 7^ fa.) 

The operator L„ is singular, because equation (|C4|l has non-trivial solutions. To discover the associated solvability 
conditions, we consider an equation of the form LnUn ~ Fn with a forcing vector Fn ~ [a-n bn Cn dn]^- The four components 
of this equation can be combined into 



{(9o' + n)[d^ + 2(2 - q)] + kld^}hn 



Akxdo{doa„ -I- 2bn) + [do + 2(2 - g)](c„ -f dod„). 



(C12) 



For the equation to be solvable, the right-hand side should not contain any term proportional to e "^^o, where uj is any root 
of the dispersion relation ((58}. If a„, b„, c„ and d„ are simply proportional to e~"^^° then the solvability condition is 

-k^uj{-iujar, + 2b„) + [-uj'^ + 2(2 - q)] (c„ - ia;d„) = 0. (C13) 

In order to satisfy the exact conditions for parametric resonance, we should consider two waves whose frequencies combine 
(by addition or subtraction) to give the frequency of the warp, which is equal to 1 in our units. Without loss of generality, 
we call the frequencies of the two waves uo and uj -\- 1. Their vertical mode numbers should be neighbouring integers, m and 
m so that they are coupled by the terms described above. The exact conditions for parametric resonance are therefore 

{-u? + m) [-o;^ 2(2 - g)] - fc^o;^ = (C14) 

and 

{-{uj + lf +m+ 1)[-{lo + 1)' + 2(2 - q)] - kl{uj + 1)^ = 0, (C15) 
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which can be satisfied simultaneously at discrete values of kx, as shown in Fig. [T] 
Thus we consider a solution in which 



(0) 



u: 



and 







LLrn 




(0) 
Win 


= A,n(n) 


1 I'm 





iuj(u)^ — m) 
(2 - - m) 



(C16) 



(0) 



-'■m + l 

"rn+1 
(0) 

^(0) 



i(a; + l)((a; + l)' 

{m+l)kx(uj + l) 
ikxioj + 1)^ 



(m - 



-I))' 
1)) 



-i(i>j + l)To 



(C17) 



while Un'^ = for other values of n, where Am{T\) and Am+i{Ti) are the slowly varying amplitudes of the two waves. 
The solution (|C16[) " (|C17I) satisfies equation (jC4[) because it is a linear combination of eigenmodes. When substituted into 
equation (|C7|) it generates components of the forcing vector J*'^^' with various frequencies. Many of these terms produce a 



non-resonant response Uli'^ that could be found by solving equation (|C7|) . Of interest here, however, are the forcing terms 



that resonate with the free wave modes. The part of Fm that is proportional to e and therefore resonates with mode m 
is (omitting viscous terms) 











(0) 




di 


Vm 

(0) 
Win 


+ 









-\kxS{m + l)v':°^ 



m+l 

\kxSmw'-^\^ 
\kxS(m + l)h^^l, 



(0) ■ 
m+l 



(C18) 



while the part of J^^'^i that is proportional to e '('^+^)^o and therefore resonates with mode m is (again omitting viscous 
terms) 



,,(0) ■ 

"m + l 

'^m+1 

(0) 
^m + 1 

."m + l. 



+ 



2 '^x ^ <^m 

1- (0) 
jlMm 



(C19) 



Applying the solvability condition HC12|) to each of these forcing vectors leads to the following relations between the amplitudes 
of the two waves: 



9lAm = 

where 
Ci 

C2 



■ClAn 



dlAn 



C2Ar. 



[uj^ioj + 2){lo + g - 1) - (2 - q){2uj + l)m](m + l)fc^ 
2(2g - 3)u}[-u}^ + 2(2 - q)m] ' 

-(2a; + l)[uj'^(uj + if - 2(2 - q)m]kl - {2q - 3)u}{u} + 2){uj'^ - m){u}'^ + 2a; + 2g - 3) 
4(2q - 3)(a; + l)[-(u; + 1)* + 2(2 - q){m + l)]fc, ' 



(C20) 

(C21) 
(C22) 



The solutions of these coupled equations are proportional to e''^'^^, with scaled growth rates si given by Si = C1C2. The true 
growth rate of the parametric instability is s = IV'i'Si. 

We have evaluated the growth rates for all the possible resonances shown in Fig. [T] (right panel). Those that involve 
couplings between an inertial mode and and acoustic mode have C1C2 < and do not lead to instability. Those that involve 
couplings between two inertial modes have C1C2 > and do lead to instability. The growth rates for q — 1.6 and = 0.01 
or 0.02 are plotted as points in Fig. [2] together with the results of the full numerical calculations. 

In the limit m S> 1 the following approximations hold: 



(15 



64(2 - g) 



m ^ + 0{m ^), 



(15 - 8g)m + 



15 -8g 



+ 0(m-i), 



C1C2 



(3-g)^(15-8g) 
1024(2 -g)2(2g- 3)- 



+ 0{m' 



(C23) 
(C24) 
(C25) 



Instability is therefore possible in this limit if g < For q = 1.6, the limiting value of the growth rate is 0.8111|7/;1. Note 
that the numerical growth rates at large kx in Fig. [2] exceed this value because of the overlap of resonances. 
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The calculation given above is readily extended to allow for a slight detuning of the resonance, of order and a slight 
damping of the waves, also of order \tp\. This leads to the system of equations 

dlAm = ClAm+l — CsAm, dlAm+1 = C2j4m — CiAm+1 + '^'^lAm+l, (C26) 

where the viscous damping coefficients are given by 

Cahw" + 2(2 - q)m] = ai{w^(-w^ + 2)kl + + l)m + 2(2 - g)[(-w^ + m)kl + m{m - 1)]} 

+i(abi - |ai){a;^[(-w^ + m)kl - m{m - 1)] + 2(2 - q)m{m - 1)} (C27) 

and a similar expression for C4, but with m 1— > m + 1 and 1— > a; + 1. Here the detuning (i.e. the mismatch between the 
frequency difference of the two waves and the frequency 1 of the warp) is wil'!/'!. In this case instability occurs at the centre 
of the resonance (wi = 0) if C1C2 > C3C4. The half- width of the resonance (in terms of frequency) can be seen from the 
condition for instability, 

2 (C3 + C4f{ClC2 - C3C4) ,„ . 
■ (<^28) 

However, in the inviscid case the width is given by 

ojf < AC1C2, (C29) 

which means that the half-width of the resonance is twice the growth rate at the centre of the resonance. In fact, the inviscid 
growth rate anywhere within the resonant width is (C1C2 — \loIY^'^ . The viscous and inviscid expressions for the width 
agree if C3 = C4 <C (CiC2)^''^. The half-width of the resonance in terms of wavenumber can be found by differentiating the 
dispersion relation: small departures and ki\il}\ from the centre of the resonance are related by 

wi ^ a;^fc^ {uj + l)^fc^ , , 

fei + 2(2 - g)m] [-(a; + 1)4 + 2(2 - g)(m + 1)] ' ^ > 

For the first successful resonance in the case q = 1.6, which has rn = 1, kx = 1.6038 and lj = —0.4374, we have 
Ci = 1.7125, C2 = 0.2692, C3 = 3.5485ai +0.2607(abi - f ai), C4 = 4.5787ai +0.7792(abi - f ai) and wi/fci = -0.3662. We 
then have instability for a < 0.2044|i/)| in the case ah = fa, or a < 0.2226|V'| in the case ab = 0. The wavenumber-half-width 
of the resonance is 3.708|'!/'| in the inviscid case. 



C2 Keplerian case 



In the Keplerian case g = | we must include viscosity in order to find laminar flow solutions. The amplitude of the laminar 
flows, and therefore the growth rate of the parametric instability, scale with IV'l/a, while the viscous damping rate scales with 
a. In order to allow these to compete, we let a = [tpl^^'^ai and ab = \ip\^^'^ahi in this section, where ai and abi are of order 
unity in the limit \tp\ <C 1. 



Now the laminar flow solutions have the expansion 



with 
C 



,1/2 



CCOST + 0{\ipf), 



,1/2 



Ssinr + 0{\i;\^), «; = 0(|V|^^^) 



1 = 0(|^|3/^), 



1 

2c^' 



S = 



1 

4qi ' 



(C31) 
(C32) 



Their phase relationship to the warp is different from the inviscid non-Keplerian case because now the epicyclic oscillator is 
driven at its natural frequency and the response is limited by viscous damping. The linearized solutions have the expansion 

Mr) = (ro, ri , • • • ) + IV^l'^'wi'' (ro, n , • • • ) + OdV-l), (C33) 

etc., where now (ro,Ti, • • • ) = (r, |V'|^^^t, • • • ). The rest of the argument goes through similarly to the non-Keplerian case, but 
with the simpler forcing vector 



-diu: 



(0) 



-ikxC cos r + (n + 1)m^+i] - C cos r w^+i + X„ 
-ifc.Ccosr[i;i°2i + (n + l)v^°li] - S sinr wl^'l, + Y„ 
-ifcxCcosr[w^°2i + + Z„ 

-ifc, C cos r[/ii°2i + {n+ l)h'^°L] 



(C34) 



The part of that is proportional to e and therefore resonates with mode m is (omitting viscous terms) 

,(0) _ Ln.„(o) 



-di 







(0) 

(0) 
Win 


+ 







-iifc.C(m + l)K^'Vi-|C7«;L+i 
-likxC{m + l)v^:^l, + liSw^°l, 
-iifeCmwf^^i 
-^ik^C{m + l)h^°l, 



(C35) 
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while the part of -F^+i that is proportional to e '('*'+i)'^o and therefore resonates with mode m is (again omitting viscous 
terms) 



(C36) 



"m + 1 






(0) 


+ 


2 in^a; ^ '-'m 


(0) 




."■m+l. 




2 Chm 



Prom the solvability conditions we now find 
(2w + + 2) - m](m + l)iA;a, 



C: 



8aiw(— w"* + m) 

{2uj + + 1)^ - m]ik:c 



(C37) 

8ai(a; + l)[-(a; + l)4 + (m + l)]- ^^^^^ 

Although Ci and C2 are now imaginary, the previous analysis can still be applied and instability is possible if C1C2 > 0, 
which again occurs for couplings between two incrtial modes. 

For the first successful resonance, which hcis m = 1, kx = 1.8795 and w = —0.4325, we have Ci = 0.1712i/Qi, C2 = 
-0.02769i/ai, C3 = 4.3750q;i + 0.2784(Q;bi - fai), C4 = 5.4174ai + 0.8608(abi - fai) and wi/fci = -0.3387. We then have 



instability for a < 0.1189|i/'|i^^ in the case Ob = |a, or a < 0.1236|i 
the resonance is 0.4065|-(/'|/q in the limit of negligible damping. 



|i/^ in the case ab = O. The wavenumber-half- width of 
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